home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / SHOOT.DEM < prev    next >
Text File  |  1991-04-29  |  3KB  |  133 lines
  1. PROGRAM d16r1(input,output);
  2. (* driver for routine SHOOT *)
  3. (* Solves for eigenvalues of spheroidal harmonics. Both
  4. prolate and oblate case are handled simultaneously, leading
  5. to six first-order equations. Unknown to shoot, it is
  6. actually two independent sets of three coupled equations,
  7. one set with c^2 positive and the other with c^2 negative.  *)
  8. LABEL 1,2;
  9. CONST
  10.    nvar=6;
  11.    n2=2;
  12.    np=2;
  13.    delta=1.0e-3;
  14.    eps=1.0e-6;
  15.    dx=1.0e-4;
  16. TYPE
  17.    gl2array = ARRAY [1..2] OF real;
  18.    gln2array = gl2array;
  19.    gl6array = ARRAY [1..6] OF real;
  20.    glnarray = gl6array;
  21.    glnvar = gl6array;
  22.    glarray = gl6array;
  23.    glindx = ARRAY [1..6] OF integer;
  24.    glinvar = glindx;
  25.    gln2byn2 = ARRAY [1..n2,1..n2] OF real;
  26.    glnpbynp = gln2byn2;
  27. VAR
  28.    c2,factr,h1,hmin,q1,x1,x2 : real;
  29.    i,m,n : integer;
  30.    delv,v : gl2array;
  31.    dv,f : glnarray;
  32.    kmax,kount : integer;
  33.    dxsav : real;
  34.    xp : ARRAY [1..200] OF real;
  35.    yp : ARRAY [1..10,1..200] OF real;
  36.  
  37. PROCEDURE load(x1: real; v: gl2array; VAR y: gl6array);
  38. (* Programs using routine LOAD must define the types
  39. TYPE
  40.    gl2array = ARRAY [1..2] OF real;
  41.    gl6array = ARRAY [1..6] OF real;
  42. and must define the global variables
  43. VAR
  44.    c2,factr : real;
  45.    m : integer;
  46. in the main routine. *)
  47. BEGIN
  48.    y[3] := v[1];
  49.    y[2] := -(y[3]-c2)*factr/2.0/(m+1.0);
  50.    y[1] := factr+y[2]*dx;
  51.    y[6] := v[2];
  52.    y[5] := -(y[6]+c2)*factr/2.0/(m+1.0);
  53.    y[4] := factr+y[5]*dx
  54. END;
  55.  
  56. PROCEDURE score(x2: real; y: gl6array; VAR f: gl6array);
  57. (* Programs using routine SCORE must define the types
  58. TYPE
  59.    gl6array = ARRAY [1..6] OF real;
  60. and must define the global variables
  61. VAR
  62.    m,n : integer;
  63. in the main routine. *)
  64. BEGIN
  65.    IF (((n-m) MOD 2) = 0) THEN BEGIN
  66.       f[1] := y[2];
  67.       f[2] := y[5]
  68.    END ELSE BEGIN
  69.       f[1] := y[1];
  70.       f[2] := y[4]
  71.    END
  72. END;
  73.  
  74. PROCEDURE derivs(x: real; y: gl6array; VAR dydx: gl6array);
  75. (* Programs using routine DERIVS must define the type
  76. TYPE
  77.    gl6array = ARRAY [1..6] OF real;
  78. and must define the global variables
  79. VAR
  80.    c2 : real;
  81.    m : integer;
  82. in the main routine. *)
  83. BEGIN
  84.    dydx[1] := y[2];
  85.    dydx[3] := 0.0;
  86.    dydx[2] := (2.0*x*(m+1.0)*y[2]-(y[3]-c2*x*x)*y[1])/(1.0-x*x);
  87.    dydx[4] := y[5];
  88.    dydx[6] := 0.0;
  89.    dydx[5] := (2.0*x*(m+1.0)*y[5]-(y[6]+c2*x*x)*y[4])/(1.0-x*x)
  90. END;
  91.  
  92. (*$I MODFILE.PAS *)
  93. (*$I LUBKSB.PAS *)
  94.  
  95. (*$I LUDCMP.PAS *)
  96.  
  97. (*$I RK4.PAS *)
  98.  
  99. (*$I RKQC.PAS *)
  100.  
  101. (*$I ODEINT.PAS *)
  102.  
  103. (*$I SHOOT.PAS *)
  104.  
  105. BEGIN
  106. 1:   write('Input M,N,C-Squared:  ');
  107.    readln(m,n,c2);
  108.    IF ((n < m) OR (m < 0)) THEN GOTO 1;
  109.    factr := 1.0;
  110.    IF  (m <> 0)  THEN BEGIN
  111.       q1 := n;
  112.       FOR i := 1 to m DO BEGIN
  113.          factr := -0.5*factr*(n+i)*(q1/i);
  114.          q1 := q1-1.0
  115.       END
  116.    END;
  117.    v[1] := n*(n+1)-m*(m+1)+c2/2.0;
  118.    v[2] := n*(n+1)-m*(m+1)-c2/2.0;
  119.    delv[1] := delta*v[1];
  120.    delv[2] := delv[1];
  121.    h1 := 0.1;
  122.    hmin := 0.0;
  123.    x1 := -1.0+dx;
  124.    x2 := 0.0;
  125.    writeln;
  126.    writeln('Prolate':17,'Oblate':23);
  127.    writeln('Mu(m,n)':11,'Error Est.':14,'Mu(m,n)':10,'Error Est.':14);
  128. 2:   shoot(nvar,v,delv,n2,x1,x2,eps,h1,hmin,f,dv);
  129.    writeln(v[1]:12:6,dv[1]:12:6,v[2]:12:6,dv[2]:12:6);
  130.    IF ((abs(dv[1]) > abs(eps*v[1])) OR 
  131.       ((dv[2]) > abs(eps*v[2]))) THEN GOTO 2
  132. END.
  133.